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一 类 种 高 频 耗 散 的 动力 学 显 式 算法 


REH, RA 


(青岛 理工 大 学 土木 工程 学 院 , 266033 青岛 ) 


摘 要 : 带 数 值 耗 散 的 算法 因 其 能 有 效 过 滤 虚 假 高 频 响应 的 影响 而 备 受 关注 。 基 于 离散 控制 理论 
2 变换 ,提出 一 类 带 数 值 耗 散 的 结构 动力 学 显 式 算 法 。 该 算法 采用 CR 法 的 速度 和 位 移 递 推 式 , 满 

足 零 振幅 衰减 率 , 且 对 线性 系统 和 非 线性 刚度 软化 系统 为 无 条 件 稳定 ,对 刚度 硬化 系统 则 是 条 件 稳 

定 的 ;该 算法 由 单个 参数 p 控制 数值 耗 散 能 力 。 通 过 对 振幅 衰减 和 周期 延长 的 理论 分 析 表 明 ,系数 

a 可 调节 算法 的 精度 和 非 线 性 稳定 区 间 ,给 出 精度 最 优 时 系数 的 取 值 。 对 特定 的 系数 取 值 该 算法 

可 转变 为 CR 法 。 通 过 算 例 对 线性 系统 和 非 线 性 系统 的 分 析 验 证 了 新 算法 具有 良好 的 精度 、 稳 定性 
和 数值 耗 散 , 表 明 新 算法 是 正确 有 效 的 。 
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A family of structural explicit algorithm with 
high-frequency dissipation 


ZHAO Jianfeng, ZHAO Xu 
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Abstract ;: Algorithms with numerical dissipation have attracted much attention because of their ability to 
effectively filter the effects of false high-frequency responses. Based on the Z-transform of discrete control 
theory ,a class of explicit time-history integral algorithm with numerical dissipation for structural dynamics 
is proposed. The speed and displacement recursion of CR algorithm are used in this new algorithm , which 
can meet the zero amplitude decay rate. It is unconditionally stable for linear systems and nonlinear stiff- 
ness softening systems ,and is conditionally stable for stiffness hardening systems. A single parameter p is 
used to control the existence of numerical dissipation and the size of the dissipation capacity. Theoretical a- 
nalysis of amplitude attenuation and period extension shows that the accuracy and nonlinear stability inter- 
val of the algorithm can be adjusted by the coefficient a , and the coefficient value is given when the accu- 
racy is optimal. And the new algorithm reduces to CR method for specific coefficients. Numerical examples 


of linear and nonlinear systems show that the new algorithm has good accuracy ,stability and numerical dis- 
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sipation , indicating that the new algorithm is correct and effective. 


Key words explicit algorithm ; numerical dissipation ; unconditional stability ; discrete control theory 


数值 积分 方法 在 求解 结构 动力 问题 上 应 用 广 
泛 , 通 过 将 运动 微分 方程 离散 为 差分 方程 求解 结构 
在 不 同时 刻 的 状态 量 ( 位 移 .速度 和 加 速度 ) 。 数 值 
积分 方法 包含 显 式 法 和 隐 式 法 两 类 。 若 当前 时 
步 的 状态 量 只 由 已 知 时 步 的 状态 量 决定 ,算法 即 为 
显 式 的 ,有 单 步 法 和 多 步 法 ;反之 则 为 隐 式 算法 。 
隐 式 算法 如 著名 的 Newmark-B iz 8 Wilson-0 法 中 
等 都 具有 很 高 的 精度 , 且 对 特定 参数 可 以 满足 无 条 
件 稳定 ,但 隐 式 算法 需要 大 量 迭 代 计 算 , 计 算 效 率 
低 ,不 适用 于 大 型 非 线 性 结构 '"""。 显 式 算法 如 中 心 
差分 法 虽然 无 需 迭 代 , 但 由 于 是 条 件 稳定 的 ,时 间 步 
长 的 选取 受到 稳定 界限 的 限制 ,过 短 的 时 间 步 长 也 
会 增加 计算 耗 时 号 ] 。 近 年 来 , 随 着 拟 动力 试验 和 实 
时 混合 试验 的 发 展 ,高 精度 、 高 稳定 性 和 高 效率 的 数 
值 积分 算法 受到 广泛 关注 ,国内 外 学 者 在 此 方面 的 
研究 取得 丰硕 成 果 。 

台湾 学 者 CHANG 提出 无 条 件 稳 定 的 显 式 算 
EU ,提高 了 显 式 算法 的 计算 效率 ,但 该 方法 仅 位 
移 表 达 式 为 显 式 。 随 后 CHEN 和 RICLESH 采用 
离散 控制 理论 分 析 NEWMARK 法 和 CHANG 方法 
后 ,推导 出 位 移 速度 双 显 式 的 无 条 件 稳定 算法 (CR 
法 )。 此 后 离散 控制 理论 在 结构 动力 学 算法 上 得 到 
PEU 。 杜 晓 琼 等 ' 国 分析 CR 法 和 CHANG 
方法 的 极点 ,提出 一 种 精度 可 控 的 无 条 件 稳 定 显 式 
算法 。GUI 等 ”总结 Newmark 算法 的 极点 ,利用 离 
散 控制 理论 推导 出 一 类 适用 于 复杂 非 线 性 问题 的 无 
条 件 稳定 显 式 算法 。TANG 等 ”| 以 全 显 式 传递 格式 
推导 出 Real-time Substructure Testing( RST) 法 ,该 方 
法 应 用 于 实时 混合 试验 仍 具 有 高 精度 |。 尽 管 以 上 
方法 都 为 显 式 且 无 条 件 稳定 ,但 不 具备 数值 耗 散 ,不 
能 过 滤 高 阶 振 型 对 结果 的 影响 。 随 后 TANG 4&7" 5| 
和 参数 op 改良 了 RST 算法 ,使 其 带 可 控 的 数值 耗 
散 ,然而 参数 的 引入 改变 了 算法 的 稳定 性 。KOLAY 
等 ("1 基于 广义 -a 法 和 CR 法 ,利用 离散 控制 理论 提 
出 无 条 件 稳 定 的 显 式 Kolay and Ricles( KR) 法 ,该 算 
法 可 以 抑制 高 阶 虚假 振 型 ,而 对 低 阶 振 型 无 耗 散 , 但 
求解 过 程 复杂 。 此 外 LI 等 ' 引 设计 了 由 单个 参数 控 
制 的 无 条 件 稳定 算法 ,具有 可 控 的 数值 耗 散 ,但 该 算 
法 为 隐 式 的 , 且 每 一 时 步 需 要 分 多 个 子 步 计算 。 

基于 离散 控制 理论 提出 一 类 具有 可 控 数 值 耗 散 


的 无 条 件 稳定 显 式 新 算法 。 本 研究 以 算法 高 频 耗 散 
为 条 件 ,利用 离散 传递 函数 推导 算法 的 参数 ,得 到 算 
法 极点 。 精 度 和 稳定 性 分 析 表 明 ,新 算法 对 线性 问 
题 和 非 线 性 刚度 软化 问题 是 无 条 件 稳定 的 ,并 且 可 
使 算法 满足 零 振 幅 衰减 率 和 极 小 的 周期 延长 率 。 最 
后 用 几 个 算 例 对 新 算法 的 性 能 进行 数值 验证 。 


1 显 式 新 算法 设计 


1.1 设计 思路 


以 线性 单 自由 度 系统 为 例 进行 推导 ,在 i 时 刻 
单 自由 度 系 统 的 运动 方程 为 
m x, + cx, + kx, = f, (1) 
其 中 :m,c, 分 别 表示 质量 .阻尼 和 刚度 ; X, ,x 
表示 系统 在 i 时 刻 的 加 速度 、 速 度 和 位 移 ; f; 表示 i 
时 刻 系 统 受到 的 外 部 激励 。 以 式 (2) 作 为 新 算法 的 
位 移 、 速 度 递 推 式 ,由 于 i+1 时 刻 的 位 移 和 速度 只 
与 i 时 刻 的 状态 量 相关 , 故 算法 是 显 式 的 。 


[^ £z Xi 十 At * X; t a AC * X; 


VM "n (2) 
Xa = X; + BAË tox, 
其 中 : a 和 8B 为 待 求 参 数 ; At 为 时 间 步 长 ; xi 和 %i, 
分 别 表示 系统 在 i+1 时 刻 的 位 移 和 速度 。 将 式 (2) 
代入 式 (1) 即 可 得 到 系统 在 i+1 时 刻 的 加 速度 。 
目前 已 有 众多 学 者 研究 离散 控制 理论 在 结构 动 
力学 算法 设计 上 的 应 用 ,但 都 是 在 已 知 算法 极点 的 
情况 下 ,进行 参数 的 推导 。 本 研究 从 不 同 的 角度 出 
发 ,以 算法 的 数值 耗 散 为 条 件 ,推导 算法 参数 ,最 终 
确定 算法 的 极点 ,完成 新 算法 的 设计 。 


图 1 新 算法 设计 思路 


Fig.1 Design ideas of new algorithm 
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L2 参数 推导 


利用 离散 控制 理论 对 式 (1)、(2) 进 行 Z 7E 
换 ” ,可 以 将 算法 的 状态 量 转化 为 关于 Z 的 二 次 方 
程 。 离 散 控 制 理论 求 得 的 系统 极点 与 算法 放大 矩阵 
的 特征 值 等 价 。 在 分 析 算 法 性 能 时 ,以 代数 运算 代 
蔡 放 大 和 矩阵 的 分 析 ,达到 简化 计算 的 效果 。 

以 系统 位 移 和 外 和 荷载 为 例 , 对 Z 变换 进行 说 明 。 
位 移 x 和 外 荷载 f 在 i 时 刻 的 Z 变换 可 表示 为 : 
Z(x) = X(z) ,2Z(f:) = F(z) ,根据 位 移 定理 “ ， 
i+1 时 刻 系统 的 Z 变换 可 由 i 时 刻 表示 , 即 Z(%;,,) = 
zX(z) 。 同 理 可 得 到 系统 速度 和 加 速度 不 同时 刻 的 
Z 变换 。 

离散 控制 理论 中 开 环 系统 的 传递 函数 由 C(z) 
表示 ,如 式 (3) ,分 子 为 系统 输出 (位 移 *) 的 Z 变换 ， 
分 母 为 系统 输入 (外 人 荷载 万 的 和 变换 。 

X(z nz + niz + no 
~ FG ] dz + diz d, Ny 
通过 对 式 (1)、(2) 的 Z 变换 ,得 到 开 环 传递 函 
数 的 各 项 系数 ,如 表 1 所 示 。 表 中 二 为 阻尼 比 ,2 = 
wAt , w 为 系统 自 振 频 率 。 系 统 特征 方程 为 

dz +diz +d =0 (4) 

SEE Zr Fe B fi. BI A Pe 3 PR RAII ei, FARE 
阵 的 特征 值 对 应 。 将 表 (1) 中 的 系数 代入 式 (4) , 解 
得 极点 表达 式 ,如 式 (5) 所 示 。 

表 1 开 环 系统 传递 函数 的 系数 


Tab. 1 Coefficients of transfer function of open-loop system 


G(z) 


n 系数 n fü d 系数 d f& 
nz 0 dz n 
ni aA? d, m(a(? +2BEQ -2) 
no (B-a) A? do ml (B-a) È -2BEQ+1] 
pol -pe0 - 92 m -4aB£Q - AB Ë -4B 
(5) 
令 极点 实 部 和 庶 部 分 别 为 n=1 - geo -È 


e= AB co -4oBRQ-4B 忆 ,那么 极点 可 表示 
为 式 (6) 的 形式 , 式 中 i 为 虚数 单位 。 
Pi =T Łe'i (6) 
算法 的 数值 耗 散 特 性 由 谱 半径 p 控制 ,p < 1 风 
表示 算法 存在 数值 阻尼 , 即 有 数值 耗 散 。 由 式 (6) 
可 得 算法 谱 半 径 为 
p= te (7) 
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由 式 (5) ~ (7) 可 知 , o 为 关于 算法 参数 w fll B 
的 函数 ,将 o 和 & 上 代入 式 (7) ,得 
p= V(B-a) -28£0 +1 (8) 
数值 阻尼 应 能 滤 除 高 阶 振 型 对 系统 响应 的 影 
Te] ,同时 对 低 阶 振 型 的 影响 应 尽量 小 ,以 保证 精度 。 
于 是 谱 半 径 应 满足 以 下 条 件 
2 一 0 时 ,p =1 (9a) 
Q— ov 时 ,p = pa (9b) 
Po 表示 2 趋 于 无 穷 大 时 谱 半 径 值 , 取 值 范围 为 
[0,1]。 显 然 , 式 (9a) 是 满足 的 。 当 2 一 x 时 , 取 
< = 0, W] 


(B-a)f +1 =p (10) 
通过 式 (10) 可 得 到 两 个 参数 之 间 的 关系 , 即 
pav! 
Pe 3 EA Sy 


式 中 ,4 为 0 的 一 次 函数 ,为 引入 结构 特性 , 取 4 = 
b(EQ 1) ,系数 /为 大 于 0 的 常数 。 参 数 a 和 8p 的 
分 母 为 (7 +4 ,将 式 (11) 代 入 式 (8) 和 式 (9a) ,得 到 
参数 8 的 表达 式 ,再 将 B 代 入 式 (11) 则 得 出 参数 a 的 
dest. 


E a*l-p 
* = P b(EQ x1) (iaa) 
P = FEOL Sa 


AIF ab 为 系数 ,将 在 第 2 节 对 其 进行 分 析 。 对 于 
参数 p。 的 选取 则 需要 考虑 系统 响应 , 若 响应 中 包含 
不 感 兴趣 的 或 虚假 的 高 频 分 量 , 则 建议 选择 较 小 的 
p-。 值 ,以 过 滤 高 频 分 量 ; 若 响应 由 低频 分 量 控制 , 则 
建议 选择 较 大 的 ps 值 ,其 取 值 应 在 [0,1 |]。 


2 算法 稳定 性 和 精度 分 析 


2.1 线性 系统 稳定 性 


线性 系统 稳定 性 通过 谱 半径 进行 分 析 , 当 谱 半 

径 p < 1 时 ,算法 是 无 条 件 稳定 的 。 式 (8) 已 给 出 新 
算法 的 谱 半 径 表 达 式 ,后 文 的 分 析 表 明 , 当 系数 a = 
b 时 ,算法 精度 最 优 。 此 时 算法 参数 如 式 (13 ) ,需要 
指出 的 是 当 a = 4 ,p。= 工时 ,新 算法 即 为 CR 法 。 
a+l - pe 


+aen+1) 


(13a) 


P 70d ee Hp) 
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注意 到 极点 虚 部 根 号 部 分 应 为 非 负 ,结合 式 
(6) 与 式 (7) 可 得 到 算法 无 条 件 稳定 时 系数 a 的 取 
值 范围 为 [1 tp. —p« V3 *p. 1 tp. tp. ' 
V3 *pi]« 
Sa = 1 将 式 (13) 代 入 式 (8) 即 可 得 到 p 与 0 
的 关系 。 如 图 2 所 示 , 阻 尼 的 存在 对 算法 的 稳定 性 
影响 ,上 且 随 着 2 的 增 大 , p 逐渐 向 p。 3E V4 OE 
够 大 时 , p = p。, 说 明 新 算法 对 线性 系统 为 无 条 件 
稳定 , 且 具 有 很 好 的 高 频 耗 散 特性 。 
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图 2 新 算法 谱 半径 


Fig.2 Spectral radius for the new algorithm 


2.2 非 线性 系统 稳定 性 


离散 控制 理论 的 闭环 传递 函数 可 以 分 析 非 线性 
系统 的 稳定 性 。 非 线性 单 自由 度 系统 的 运动 方程 改 
为 增 量 形式 如 式 (14) , 式 中 AL; = Af, -Axi,k, 为 
非 线 性 系统 i 时刻 的 刚度 ,随时 间 不 断 变 化 。 


mA X; + CA X; = AL; (14) 
闭环 系统 的 传递 函数 可 表示 为 
e s AX(z) |  G'(z) (15) 


X AF(z) 1+kG’(z) 
其 中 G'(z) 表达 式 为 


ta ; ; 
nz 十 nz * ng 


RE dz * d,z +d, (16) 
式 (16) 中 的 各 项 系数 见 表 2。 
于 是 闭环 系统 特征 方程 为 
ELLA (17) 


2+dz+d 
R2 闭环 系统 传递 函数 的 系数 


Tab.2 Coefficients of transfer function of close-loop system 


n 系数 n fü d 系数 d fü 
n; 0 d, m 
n, BAP d, 2m(B&Q -1) 
no (ay -B)AP do m(1 一 2650) 
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工程 中 常见 的 非 线性 问题 包括 刚度 软化 
(0 « k, <k) FILII BE i f (k, >k) 两 种 ,定义 如 图 3 
所 示 。 


图 3 刚度 非 线性 定义 
Fig.3 Definition of stiffening and softening behavior 
利用 MATLAB 绘制 非 线 性 系统 的 根 轨迹 图 , 当 
轨迹 在 单位 圆 上 或 以 内 则 系统 是 稳定 的 。 取 a = 
b=1, p, =1,At=0.2s,w = 7,m=1,é =0, 绘 
制 根 轨迹 如 图 4 所 示 。 
LS 


-20 -1.5 -10 -05 0 05. q10- 15 
实 轴 


图 4 GC) 根 轨迹 图 
Fig.4 Root locus for G(z) 

图 4 中 共有 两 条 根 轨迹 ,一 条 起 于 p =0 经 pr 
止 于 zi , 男 一 条 起 于 p, 2p IET z = % pl， 为 
k, = 0 时 特征 方程 的 极点 pia W k, = 天 时 特征 方程 
的 极点 。 刚 度 软化 系统 对 应 的 根 轨迹 为 pl pis, 
全 部 在 单位 圆 上 ;刚度 硬化 系统 对 应 的 根 轨迹 为 
Pia 一 212 ,在 z =- 1 处 穿 出 单位 圆 。 故 算法 对 于 
非 线 性 刚度 软化 系统 为 无 条 件 稳 定 , 刚 度 硬 化 系统 
为 条 件 稳定 。 

将 z = -1 RARO?) ,可 得 刚度 硬化 系统 的 稳 
定性 区 间 如 式 (18) ,可 以 看 出 稳定 界限 与 系数 a 的 
取 值 成 反比 而 与 p。 的 取 值 成 正比 。 对 于 无 阻尼 系 
统 , 非 线性 稳定 区 间 只 与 a 有关 。 
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2.3 算法 精度 分 析 


一 个 数值 算法 的 精度 可 由 振幅 衰减 率 (AD ) 和 
周期 延长 率 ( PE) 两 个 指标 衡量 ,定义 如 下 。 
AD=1 - exp( - 2m £) (19) 
PE-(0/0-1 (20) 
其 中 : É = -In(g? t2)/0,0 = arctan( e/o.) 。 由 图 
2 可 知 , =0,p。= 1 时 , 谱 半 径 恒 为 1, 即 o +e” = 
1, 此 时 上 =0 ,也 就 是 AD =0 。 因 此 以 8=0,p。 =1 
分 析 新 算法 的 PE 并 探讨 系数 a 的 取 值 。 
首先 分 析 系 数 取 任意 值 时 新 算法 的 PE, 考虑 
a=1,b=1~4 和 a = 2,b = 1 ~4 两 种 情况 , 绘 
H PE 图。 如 图 5 所 示 , 两 种 情况 下 仅 当 a = 时 PE 的 
绝对 值 最 小 , 当 ae b 时 ,尽管 A/T 取 值 足够 小 ,PE 
绝对 值 仍 然 很 大 .图 6 ~ 7 为 a=5 时 的 PE 图 ,图 6 可 
以 看 出 ,周期 延长 率 的 绝对 值 随 系 数 的 增加 呈现 先 增 
后 减 的 变化 规律 ,在 9 ~ 13 范围 内 精度 较 高 .图 7 为 
a =b =9 ~ 14 的 范围 内 取 值 的 PE 图 ,在 a=6b=11 时 ， 
周期 延长 率 的 绝对 值 最 小 ,不 足 0. 002 ,精度 最 高 。 
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5 任意 系数 的 PE 


Fig.5 PE values corresponding to arbitrary coefficient 
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图 6 a = b W PE(1) 
Fig.6 PE values corresponding to a = b (part 1) 
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图 7 a = 时 算法 PE(2) 
Fig.7 PE values corresponding to a = b (part 2) 
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图 8 =b = 11 时 算法 AD 
Fig.8 AD values corresponding to a = b = 11 
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图 9 不 同 阻尼 比 的 算法 AD 
Fig.9 AD values with different damping ratio 


结合 以 上 分 析 , 取 a =5=11 绘制 新 算法 的 AD 
图 ,由 图 8 可 知 ,新 算法 AD s =1 时 最 小 且 恒 为 
0 ,这 与 前 述 分 析 一 致 。 新 算法 的 AD 与 p。 呈 负 相 关 ， 
ps 的 取 值 越 接 近 1,AD 的 变化 越 明 显 。 图 9 进一步 展 
现 阻尼 比 对 算法 AD 的 影响 ,上 =0 时 AD ENG, pE 
着 阻尼 比 的 增加 AD 逐渐 增 大 。 因 此 , 系数 取 值 为 
a=b=11 上 且 p。=1 时 ,算法 精度 最 高 ,但 此 时 新 算 
法 无 数值 阻尼 。 最 终 参 数 a MB 中 仅 含 fI pu 两 个 
常数 系数 , 即 式 (13 ) 。 
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3 多 自由 度 系统 参数 


多 自由 度 系 统 在 i 时 刻 的 运动 方程 为 
Mx,+Cx,+Kx,=F, (21) 
式 中 :MC 和 分 别 为 系统 的 质量 .阻尼 和 刚度 和 矩 
阵 ; xx, RI x; 分 别 为 系 加 速度 .速度 和 位 移 向 量 ; 
F, 为 外 荷载。 
多 自由 度 系统 算法 参数 与 单 自 由 度 系统 的 式 (2) 
和 式 (12 ) 相 似 , 带 可 挖 数值 阻 尼 新 算法 的 递 推 式 为 


Xi = X, At* x, +A - X, 
| MA (2) 
Xi, = Xj +BAr ”Xi 
其 中 参数 为 
a - (a*1-p)-*(a-M *a/2At-C € A -K) M 
(23) 


B=a: (a-: M+a/2At - C« AC -K) ! -M (24) 


4 算 例 分 析 


本 节 通 过 2 个 算 例 验证 所 提 新 算法 对 线性 和 非 
线性 结构 动力 分 析 的 精度 、 稳 定性 和 数值 耗 散 性 能 。 

算 例 1 单 自 由 度 弹 得 质量 阻尼 系统 ,其 质量 
为 m =1 ,刚度 为 k =4m” ,分 别 取 上 = 0, £ 20.02, 
初始 条 件 x(0) = 1.0 m,v(0) = 0。 分 析 位 移 响应 。 
选 CR 法 和 中 心 差分 法 (CD) 与 新 算法 ( 取 a = 11, 
pa ,图 中 为 Designed) 对 比 , 取 Al = 0.1s ,图 中 
Analytical 为 解析 解 。 

图 10 为 系统 的 位 移 响 应 ,新 算法 不 存在 周期 误 
差 累 计 , 周 期 延长 率 在 -4.3% 左右 ,与 解析 解 吻合 
较 好 。 
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Analytical =-=-Designed +- CR= =CD 
1.0 à q è . t 5 es 
jl 
s | 
Th 
a 
aA 0.0 
-0.5 
-1.0 
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i 时 间 /s i 
图 10 系统 位 移 时 程 响应 (上 =0) 
Fig. 10 Displacement time histories of system with £ =0 
而 CR 法 和 中 心 差 分 法 存在 周期 误差 累计 , CR 
法 为 周期 延长 ,平均 周期 延长 率 为 11. 1% , 极 值 为 
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27.0% ;中 心 差 分 法 则 为 周期 缩短 ,平均 周期 延长 率 
H -8.6% , 极 值 为 -17.0% ,上 且 随 着 时 间 增 长 ,周期 
误差 仍 有 增 大 的 趋势 。 图 11 中 规律 与 图 10 相同 ， 
阻尼 的 存在 对 新 算法 无 周期 累计 误差 影响 ,周期 延 
长 率 在 -3.9% 左右。 


ss 


Analytical 一 -一 -Designed =- CR= -CD 
1.0 F 1 È i E 3 
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E 
bo 
ij 00 
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1 时 间 /s i 
图 11 系统 位 移 时 程 响 应 (和 =0.02) 
Fig.11 Displacement time histories of system with é = 0. 02 
算 例 2 分 析 新 算法 的 高 频 耗 散 特 性 。 双 自由 
度 体系 的 运动 方程 如 式 (25) 信 ,mi 2 m, = 1,h = 
10*,k, = 1 , 自 振 频 率 w， = 1.0,w，= 100.0, 各 参数 
均 为 国际 单位 ,底层 位 移 解 析 解 分 别 为 式 (26) 和 式 
(27) ,wow 考虑 高 频 响应 , Xotorn 为 不 考虑 高 频 响 
应 的 结果 ,其 对 应 的 速度 为 vs 和 vs。。 取 At = 
0.02 s ,用 本 研究 算法 不 同 p。 的 结果 与 解析 解 
对 比 。 


m, 0f% ki vk, -kq[w k,sin(1. 2t) 
| 0 Ill - k, k, hai | 0 | 
(25) 

Xaa 7 8in(1. 2t) +0. 000 3sint —0. 012sin( 1007) 
(26) 
X elernc = Sin(l. 24) + 0. 000 3sint (27) 
图 12 中 新 算法 位 移 计 算 结果 与 解析 解 吻合 较 
好 ,但 从 图 13 局 部 放大 图 可 以 看 出 高 频 耗 散 存 在 与 
否 的 差别 。 p。= 1.0 时 无 耗 散 , 计 算 结果 受 高 频 影 
响 明 显 ,与 xu 存在 较 大 偏差 ,而 p。= 0 时 算法 高 
频 耗 散 强 ,此 时 的 位 移 与 xau. 完全 重合 。 算 法 数 
值 耗 散 能 力 在 位 移 时 程 中 表现 得 更 明显 。 对 比 图 
14(a) 和 14(b) 的 速度 时 程 可 以 看 出 ,数值 阻尼 的 存 
在 也 使 得 新 算法 计算 结果 与 vw 重合 , 且 随 着 p. 
的 减 小 , 即 数 值 阻尼 的 增 大 ,新 算法 计算 结果 向 解析 
解 靠 近 越 快 ,表明 耗 散 能 力 越 强 。 在 图 14(¢c) 中 表现 
得 尤为 明显 , p。= 0.8 时 ,新 算法 耗 时 1. 5 s 才 消除 
高 频 响应 的 影响 ,而 p。= 0 时 , 仅 需 0.6 s。 由 此 可 
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见 ,新 算法 有 良好 的 高 频 耗 散 能 力 。 
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Fig.12 Time histories of displacement response 
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Fig. 13 Local zooming view for time histories of displacement 
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Fig.14 Time histories of velocity response 
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基于 离散 控制 理论 ,采用 CR 法 的 位 移 和 速度 
递 推 式 ,提出 一 类 带 数 值 耗 散 的 结构 动力 学 算法 。 
结论 如 下 。 

1) 新 算法 为 显 式 算法 ,对 线性 系统 和 非 线性 刚 
度 软化 系统 是 无 条 件 稳定 的 ,而 对 非 线 性 刚度 硬化 
系统 是 条 件 稳定 。 

2) 系数 a 可 控制 算法 精度 , 随 着 系数 值 的 增加 
算法 精度 呈 先 增 后 减 的 变化 规律 ,在 a = 11 ,p。= 
1 时 精度 最 高 ,可 满足 零 振幅 衰减 率 ,周期 延长 率 小 
于 0.002。 

3) p, 控制 算法 的 数值 阻尼 , o. = 1 时 算法 无 
数值 阻尼 ,p。= 0 时 数值 阻尼 最 大 ,高 频 耗 散 最 强 。 

4) 算 例 分 析 验 证 了 新 算法 的 精度 和 稳定 性 ,对 
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单 自 由 度 系统 周期 延长 率 约 为 -3.9% , 优 于 中 心 差 
分 法 和 CR 法 。 
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